classdef Bessel < FilterFather 
    properties
        NormalizedPoles = [];
    end
    
    methods
        function this = Bessel(guiData)
%             disp('holi')
            this.calculateFilter(guiData);
            minEdge = min([guiData.Fp guiData.Fa])*2*pi;
            maxEdge = max([guiData.Fp guiData.Fa])*2*pi;
            if this.Order < 30
                this.calcTimeData([1,2]);
                this.calcData([minEdge maxEdge]);
            end
%             disp('holu')
        end
        
        function calculateFilter(obj, guiData) % Load GUI information.
            obj.GroupDelay = guiData.Delay;
            Ap = guiData.Ap;
            Fp = guiData.Fp;
            Aa = guiData.Aa;
            obj.DesnormType = guiData.DenormType;
            obj.DesnormValue = guiData.DenormVal;
            Fa = guiData.Fa;
            Tolerance = guiData.Tolerance;
            obj.FilterType = guiData.filterType;
            MaxOrder = guiData.order;   % Default is zero. LCDTMAB!!!!
            MaxQ = guiData.maxQ;
%             MaxOrder
            if MaxQ && ~MaxOrder
                MaxOrder = 9;
%                 disp('holu')
            end
%             switch (guiData.DenormType)
%                 case 'Fp'
                    fdenorm = Fp;
%                     obj.WaNormalized=Wa/Wp;
%                     obj.Wdesnorm = Wp;
%                 case 'Fa'
%                     fdenorm = Fa;
%                 otherwise
%                     fdenorm = Fp;
%             end
%             if (MaxOrder == 0)
%                 MaxOrder = 15;    % Change nonsense default values to
%             end
            if (MaxQ == 0)
                MaxQ=30;        % sense default values.
            end
            
            if (MaxQ<0.5)
                MaxQ=0.5;
            end
            n = 0;
            switch obj.FilterType
                case 'Lowpass'
                    [n,~,~,q] = obj.fitAttenuation(Fp,Ap,Fa,Aa,MaxQ,MaxOrder,fdenorm);                
                case 'Group-Delay'    % GD, Tolerance GD are specified.
                    guiData.Fa = 10*guiData.Fp;  %If not, when plotting, starting edge = -inf
                    [n,~,~,q] = obj.fitGroupDelay(Fp,obj.GroupDelay,1-Tolerance/100,MaxQ,MaxOrder);
                otherwise
                    errordlg('Bessels aproximations does not contemplate that type of filter.','Inconsistent Error')
            end
            % On error: fitAttenuation returns n=0.  
            if ((n > 0)&&(n < 10)) % 1<=n<10
                obj.Order = n;
                obj.Poles = obj.NormalizedPoles;
                obj.K = abs(prod(obj.Poles));
                obj.Zeros = [];
                obj.frecDesnorm(Aa);
                obj.K = abs(prod(obj.Poles));
                obj.desnormalization();
%                 obj.TransferFunction = tf(num,den);
                obj.TransferFunction = obj.TransferFunction*10^(guiData.gain/20);
            else
                obj.Order = 40;    %Error code.
                obj.K = 0;
                obj.Poles = 1;
                obj.Zeros = [];
            end
        end

        function pol = BesselPoly(obj,N)
        % Returns the polynomial of Nth order Bessel.
        % pol = BesselPoly(N)
        % Bn is the Sum from k=0 to N of {b(k)*s^k}
        % b(k) = (2N-k)! / k!(N-k)! 2^(N-k) ; k=0,1,..,N
            b = zeros(0,N+1);
            for k = 0:1:N
                b(k+1) = factorial(2*N-k) / (factorial(k) * factorial(N-k) * 2^(N-k));
            end
            pol = fliplr(b);   % The fliplr is to sort values in ascendent powers of S.
        end
        
        function [z,p,k] = BesselPZ(obj,N)
            p = obj.BesselPoly(N);    % Returns the Bessel poly Nth order.
            z = 1; k = 1;	% Bessel filters do not have finite zeros.
            p = p / p(end); % Normalize.
        end
        
        function w=LookForAttenuationFrecuency(obj,N,Ap,wi)
        %Gets the frecuency of attenuation=Ap in a Nth order Bessel polynomial.
        %w=LookForAttenuationFrecuency(N,Ap,wi)
            [z,p,k] = obj.BesselPZ(N);	% Returns the [z,p,k] of Bessel filter Nth order.
            H = k*tf(z,p);          % Bessel Transfer Function.

            % The numerical equation is solved by the secant method:
            % http://en.wikipedia.org/wiki/Secant_method
            Nmax=50;n=1;    % Maximum iteration number.
            w0=0;w1=min(2,wi);     % Starting frequencies. w must be in (0,wi).
            w=w1-(w0-w1)*(bode(H,w1)-10^(-Ap/20))/(bode(H,w0)-bode(H,w1));
            while ((abs(w-w1)>10e-5)&&(n<Nmax))
                w1=w0;
                w0=w;
                if ((w0>=0)&&(w1>=0))
                    w=w1-(w1-w0)*(bode(H,w1)-10^(-Ap/20))/((bode(H,w1))-(bode(H,w0)));
                    n=n+1;
                else
                    w=0;    % Error code.
                    break;
                end
            end
        end
        
        function [Ret,num,den,q]=checkAttenuation(obj,N,fp,Ap,fa,Aa,fdenorm)
        % [Ret,num,den]=checkAttenuation(N,fp,Ap,fa,Aa)
        %   |->Checks if satisfy Low Pass template.
        % Attenuations are in dB.
        % If Ret is false, only it is returned. Does not return num or den.
            p=obj.BesselPoly(N);	% Returns the polynomial of Nth order Bessel.
            obj.WaNormalized = fa/fp;
            
            wp=obj.LookForAttenuationFrecuency(N,Ap,obj.WaNormalized);
            
            if (wp==0)
                Ret=false;
                errordlg('Reached physical limitations. Bessel filter cannot fit the specifications.','Physical limitations');
                return;
            end

            %Gets the frecuency of attenuation=Ap in a Nth order low-pass Bessel tf.
            p=p/p(end);                 % To normalize the denominator.
            [num,den]=lp2lp(1,p,1/wp);  % To have normalized freqeuency wp=1 => Ap.
            obj.Poles = roots(den)';
            obj.Zeros = [];
            obj.K = 1;%abs(prod(obj.Poles));
            obj.NormalizedPoles = obj.Poles';
            den = poly(obj.Poles);
            H=tf(num,den);
            
            [MagA,~] = bode(H,fa/fp);   % Evaluate H in Wpnormalized=fa/fp.
                
            [num,den]=lp2lp(num,den,fdenorm*2*pi);    % Desnormalize
            obj.Wdesnorm = fdenorm*2*pi;
            Ret = ((-20*log10(MagA))>=Aa);

            obj.calcQ();    % Calc the Q of all poles, orderer.
            q=max(obj.Q);   % Here is QMax.
        end
        
        function [n,num,den,q]=fitAttenuation(obj,fp,Ap,fa,Aa,MaxQ,MaxOrder,fdenorm)
        %Returns n, num and den for specified parameters on Bessel aprox.
        %[n,num,den]=fitAttenuation(fp,Ap,fa,Aa)
        % |-> n, num and den are calculated given a Low Pass template.
            Nmax=10;n=1;
            [Ret,num,den,q]=obj.checkAttenuation(n,fp,Ap,fa,Aa,fdenorm);
            while ((~Ret) && (n<15) && (~MaxOrder) )
                n=n+1;
                [Ret,num,den,q]=obj.checkAttenuation(n,fp,Ap,fa,Aa,fdenorm);
            end
%             if (q>MaxQ)
%                 if (n>1)
%                     [~,num,den,q]=obj.checkAttenuation(n-1,fp,Ap,fa,Aa,fdenorm);
%                     n = n - 1;
%                 else
%                     errordlg('Minimum order reached. Filter cannot fit the specifications','Minimum order reached');
%                 end
%             end
            
            if (n>=Nmax)
                num=0;den=1;n=0;        % Error code.
                errordlg('Reached physical limitations. Bessel filter cannot fit the specifications.','Physical limitations');
            end
        end
       
        function [n,num,den,q]=fitGroupDelay(obj,fp,GD,tol,MaxQ,MaxOrder)
        %Returns n, num and den for specified parameters on Bessel aprox.
        %[n,num,den]=fitGroupDelay(fp,GD,tol,MaxQ,MaxOrder)
        % |-> n, num and den are calculated given a GD specifications.
%             disp('holu')
            Nmax = 10;
            if MaxOrder
                n = MaxOrder;
            else
                n = 1;
            end
            
            [Ret,num,den,q] = obj.CheckGroupDelay(n,fp,GD,tol);
            
            while ( ((~Ret) && (n < 15)) && ~MaxOrder )
                n = n + 1;
                [Ret,num,den,q] = obj.CheckGroupDelay(n,fp,GD,tol);
            end
            
%             if (q > MaxQ)
%                 if (n > 1)
%                     [~,num,den,q] = obj.CheckGroupDelay(n-1,fp,GD,tol);
%                     n = n-1;    
%                 else
%                     num=0;
%                     den=1;
%                     n=0;        % Error code.
%                     errordlg('Minimum order reached. Filter cannot fit the specifications','Minimum order reached');
%                 end
%             end
            
            if (n >= Nmax)
                num = 1;
                den = 1;
                n = 40;        % Error code.
                errordlg('Maximum order reached. Filter cannot fit the specifications','Maximum order reached');
            end
        end
        
        function [Ret,num,den,q]=CheckGroupDelay(obj,n,fp,tau,tol)
        % [Ret,num,den]=CheckGroupDelay(n,fp,tau,tol)
        % tau is the groupdelay given by user.
%             disp('holi')
            den = obj.BesselPoly(n);% Returns the Bessel poly Nth order.
            num = 1;                % Bessel filters do not have finite zeros.
            den = den/den(end);     % Normalize.
            H = tf(num,den);        % Bessel Transfer Function.

            waux = tau*2*pi*fp;       % Change of variables.
            Ba = factorial(2*n)/(factorial(n)*(2^n));             % Calculus of
            Delay = 1 - (((bode(H,waux))^2)*(waux^(2*n))/(Ba^2));   % delay
            [num,den] = lp2lp(num,den,1/tau);
            fdenorm = fp;

            obj.Poles = roots(den)';
            obj.Zeros = [];
            obj.K = 1;%abs(prod(obj.Poles));
            obj.NormalizedPoles = obj.Poles';
            den = poly(obj.Poles);
            
%             [num,den] = lp2lp(num,den,fdenorm*2*pi);    % Desnormalize
            obj.Wdesnorm = fdenorm*2*pi;
            Ret = (Delay>=tol);
            
%             obj.Wdesnorm = fdenorm*2*pi;       
            
            obj.calcQ();    % Calc the Q of all poles.
            q = max(obj.Q);   % Here is QMax.
        end
    end
end

